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^^ . As a first step in the first passage problem for passive tracer in 

Q>^ I stratified porous media, we consider tlie case of a two-dimensional sys- 

tem consisting of two layers with different convection velocities. Using 
a lattice generating function formalism and a variety of analytic and 
numerical techniques, we calculate the asymptotic behavior of the first 
'^ ' passage time probability distribution. We show analytically that the 

5 ' asymptotic distribution is a simple exponential in time for any choice 

Q , of the velocities. The decay constant is given in terms of the largest 

eigenvalue of an operator related to a half-space Green's function. For 

the anti-symmetric case of opposite velocities in the layers, we show 

that the decay constant for system length L crosses over from L~^ 

C^ I behavior in diffusive limit to L~^ behavior in the convective regime, 

where the crossover length L* is given in terms of the velocities. We 
also have formulated a general self-consistency relation, from which 
we have developed a recursive approach which is useful for studying 
the short time behavior. 

Keywords: First passage problem; convection-diffusion equation; layered 
system; asymptotic behavior. 



1 Introduction 

The motion of a passive tracer in a fluid under the combined action of molec- 
ular diffusion and convection arises in a variety of settings, such as fluid flows 
through porous media, fixed-bed catalytic reactors, and the dispersion of pol- 
lutant in oceans |jl|. In many situations, the convection-diffusion equation 
(CDE) describing the variation of tracer concentration with space and time 
becomes inhomogeneous, i.e., the fluid velocity field and/or the diffusivity 
is not a constant, but a function of spatial position. One obvious method 
in the study of inhomogeneous systems is a perturbation technique 0, §]. 
Here, one starts from a homogeneous version of the system, which usually is 
solvable. The velocity or diffusivity is written as a sum of a homogeneous 
and an inhomogeneous term and the appropriate quantities are expressed as 
expansions involving the inhomogeneous term. The perturbative method is 
not immediately applicable when the magnitude of the disorder is strong, 
as many or an infinite number of terms are required, but in some cases, 
effective medium theories may be used for an approximate summation [^]. 
However such techniques are certainly not suitable for systems whose disorder 
is strongly correlated in space, where usually only straightforward numerical 
simulation or, for problems with an appropriate geometry, network models 
P, are useful. 

In this paper we are motivated by the particular case of tracer dispersion 
in porous geological systems such as aquifers and hydrocarbon reservoirs P, 
0] , and by the observation that such materials are very prominently stratified 
[^]. In this context, Matheron and de Marsily [^ first observed that when 
the number of layers is effectively infinite, the velocity fluctuations associated 
with the variation in structure and permeability of the layers could give rise 
to superdiffusive tracer motion. Several authors studied this problem further 
[ Pm , |n|, |T^, and by now there is a fair understanding of the tracer probability 
distribution for the case of a large number of horizontally-infinite layers. 
Unfortunately, the results to date do not provide concrete statements about 
the most practical configuration, involving a source and sink of tracer at finite 
separation. One would like to solve the first passage time problem for a large 
number of horizontal layers of finite extent, with various boundary conditions 
(sink or reflection) at the system edges. As a first step in this direction, we 
consider the simple case of tracer motion in a geometry consisting of two 
two-dimensional, semi-infinite layers, where tracer is released in the interior 



point and is adsorbed at the edges. Although a great simphfication compared 
to the case of an infinite number of layers, as we shall see this problem is 
already sufficiently difficult that only an approximate solution is available. 
(In fact, even in the ostensibly elementary problem of simple diffusion in 
two half-spaces with different diffusivities, a lengthy analysis has recently 
appeared ||13| .) 

The analysis to follow is based on an exact generating function formal- 
ism for biased random walks in the geometry of interest, and approximation 
schemes to extract the asymptotic behavior. More generally, we hope that 
our methods are pertinent to the problem of transport in system with "block" 
disorder, for inhomogeneous materials which are naturally modeled as a col- 
lection of finite homogeneous sub- regions placed in contact WM . When the 



size of the sub- regions is much less than that of the system itself, or the wave- 
length of any probe, the disorder is short ranged and perturbation techniques 
are appropriate, but otherwise few methods beyond numerical simulation are 
available. 

In this paper, we address the first passage time properties of passive tracer 
which convects and diffuses in a two layer system, which is the simplest non- 
trivial case of layered structures. The system, shown in Fig. 1, consists of 
two semi-infinite blocks occupying the two-dimensional region |x| < L. The 
blocks are in physical contact, allowing tracer to pass between them, and 
inside each block different fluid flow flelds convect the tracer. For simplicity, 
the tracer diffusivities are assumed to be equal. The two flnite boundaries at 
X = ±L are taken to be perfect absorbers. Tracer is released at some point in 
the interior, and the time-dependent flux at the boundary is computed, which 
in this situation is identical to the flrst passage time probability distribution. 

We begin in Section 2 with a precise formulation of the model as a random 
walk process, and by introducing appropriate generating functions, formu- 
late an exact self-consistency relation for the flrst passage time distribution. 
In Section 3, in order to obtain the asymptotic behavior at long time, we 
expand the flrst passage time distribution terms of the number of times a 
tracer particle has crossed the interface between the blocks before reaching 
the boundaries. We then approximately sum the expansion, using the central 
limit theorem, to obtain the asymptotic distribution. We show analytically 
that the asymptotic distribution decays as a simple exponential in time, for 
any choice of the velocity fields, where the decay constant is given in terms of 
the largest eigenvalue of an operator which is related to a half-space Green's 



function. We estimate the decay constant for the special case of the "an- 
tisymmetric" modeL For the hmiting case of high velocities, we estimate 
the largest eigenvalue and find the decay constant behaves as 1/ L, which 
agrees with numerical simulations. In the opposite case of pure diffusion, 
the decay constant behaves as l/L^, in good agreement with analytic esti- 
mates and numerical simulations. In Section 4, we consider the behavior in 
the intermediate velocity regime, using two methods: an expansion method 
about the convective limit, and a more general scaling argument which pre- 
dicts a crossover from a diffusive to a convective regime as L increases. The 
crossover length L* is given in terms of the velocity, and the scaling argument 
is consistent with the above results as well as those of numerical simulations. 
We conclude in Section 5, with a summary and discussion of future possibil- 
ities. In Appendix A, we interpret the general self-consistency condition as 
a recursion relation, and obtain an expansion useful for obtaining the short 
time behavior of the first pasage time distribution. Appendix B gives solves 
the first passage time problem explicitly for the simple case of convection 
and diffusion in a single layer. 

2 Self-consistency Relation 

2.1 Definition of the Model 

Since the tracer motion is given by the convection-diffusion equation, one 
may equivalently think of it as a biased random walk on a spatial lattice 
in discrete time. Consider then a lattice of unit spacing in the x-y plane, 
where the velocity field takes on different values in the upper and lower half- 
planes, and where only the region — L < a; < L is relevant - see Fig. 1. The 
probability P„(x, y) that the particle is at position (x, y) at time n is given 
by the master equation 

Pn+i{.x,y) = p^{y-l) py{y -1) Pnix-l,y-l) 

+ p,{y+l)[l-py{y + l)]Pn{x-l,y+l) 
+ [l-p,{y-l)]py{y-l) Pn{x + l,y-l) 

+ [i-pAy + 1)] [1 - Pyiy + 1)] Pn{x + i,y + i) 



where Px (py) are the hopping probabihties in the positive x (y) direction, 
which satisfy 

Px/y{y) j pd^^ otherwise, ^^ 

and where the Kroneker deltas prescribe that the particle starts from (xo, Vo) 
at time n = 0. The master equation implies that each step is along the 
diagonal of a square, which is a particularly convenient hopping rule for the 
analysis to come, and in the limit of long time and distance, as good as 
any other. Indeed, by expanding the right hand side in a Taylor series about 
(x, y), it is easy to see that (|l|) is equivalent to a convection-diffusion equation 
with diffusion coefficient 1/2 and velocity {2px — 1, 2py — 1). (There are also 
higher order terms involving the derivatives of px/y which are not relevant in 
the cases considered subsequently). 

The first passage problem corresponds to absorbing boundaries at the 
system edges, so we put Pn{x,y) = at a; = ±L, and define Hi^ to be the 
probability that the particle first reaches x = ±L at time n. Motivated by 
simplicity, and previous work on the many-layer problem, we suppose the 
velocities are in the x-direction, parallel to the layer boundaries, so that 
Py=Py = \- For the same reasons, we assume that the net convective bias 
or average velocity vanishes, which implies that the probability of hopping 
to the right in the upper half plane equals the probability of hopping to 
the left in the lower half plane, or p" = 1 — p^. We refer to this as the 
"antisymmetric" model, and some results about the general case appear in 
Section 5. Lastly in the remainder of this section, we simplify the analysis by 
further assuming p" = 1 (and pf = 0); so that particles in the upper (lower) 
half plane always move to the right (left), and we refer to this as the "+/— 
model." The latter restriction is lifted in Section 4. 

2.2 Derivation of a self-consistency relation 

In this section, we derive some useful relations for the +/— model. The 
master equation (|l|) reduces to 

r (P„(x-l,y-l) + P„(x-l,y + l))/2 y>2 
Pn+i{x,y) = {Pn{x-l,y + l) + Pn{x + l,y-l))/2 y = 0,l 
[{Pn{x + l,y-l) + Pn{x + l,y + l))/2 y<-l 

+ ^n+l,0 ^x,xo ^y,yo- (3) 



We define the following generating functions 

oo 

n=0 

oo 

G+{x,a,z) = ^P{x,y,z)ay, 

y=l 


G-{x,a,z) = Y. Pix,y,z)ay. (4) 

y=~oo 

Substituting Eq. (^ in Eq. (^), we obtain 

z 1 

G+(x,a,z) = -(a-\ — )G+(x — l,a,z) 

2 a 

+ ^{P{x + 1,0, z)a-P{x- 1,1, z)) 

+ S^,xoay° (5) 

and 

z 1 

G_(x,a,z) = -{a -\ — )G_(x + l,a,z) 

2 a 

+ ^{Pix ~ 1,1, z)-P{x + 1,0, z)a), (6) 

where we assume yo > 1 without the loss of generality. 

The functions G+ and G- can be expressed in terms of simple Green's 
functions. We define the Green's function in the upper block {y > 1) to be 
the solution of 

z 1 

g+{x; x' , a, z) = -(a + -)5f+(x - 1; x' , a, z) + 4,^/, (7) 

z a 

which is 

, , X f (f(a + ^))"""' ii-L<x'<x<L ,^, 

^+^^'^'«'^) = |o " otherwise. " ^^^ 

Similarly, the Green's function for the lower block {y < 0), the solution of 

z 1 

g_{x] x' , a,z) = -{a\ -)g-{x + 1; x', a, z) + 5^,^/, (9) 

6 



which is 

g_[^x,x,a,z} |q otherwise. ^^> 

Using these Green's functions, G± can be expressed as 

z ^-^ 
Gj^{x,a,z) = 77 X! g+{x;x\a,z){P{x' + l^Q^ z)a — P{x' — 1^1^ z)) 

^ x>=-L+l 

+ g+{x]Xo,a,z)ay°, (11) 

and 

z ^-^ 
G-{x,a,z) = - Y. 9-{x;x',a,z){P{x'-l,l,z)-P{x'+l,0,z)a) {12) 

^ x'=-L+l 

There is a simple way to understand Eqs. (0) - ([l^) . The Green's function 
g+ {gj) is the solution of the homogeneous convection equation Px = ^ (0) for 
a particle starting from (x',0). The second term in Eq. (^TJ) corresponds to 
the original particle which starts from (xo, Ho)- The first term of the equation 
is due to the existence of the boundary. It subtracts the contribution of the 
particle [P{x, 1, z) term] which leaves the block, and add the contribution of 
the particle [P(x, 0, z) term] which enters the block. The equation for the 
lower block Eq. (0) has essentially the same structure, except it lacks the 
second term due to the absence of a starting particle in the block. 

We are interested in the first passage properties which can be calculated 
from H^. They are related to G± by 

H+{z) = z G+{L- 1,1, z) 

H-{z) = z G4-L + 1,1, z), (13) 

where H^{z) is defined to be E^=o^n^"- In Eqs. (|llD - (|12D, G± are ex- 
pressed in terms of two unknown functions P{x,0,z) and P{x,l,z), which 
again can be calculated from G± themselves as follows. We expand Eq. (p!l|) 
as a series of a. The terms proportional to a, from the definition of G+(x, a, z), 
are exactly P{x, 1, z) a. In other words. 



X I I 



x-Xo + Vo- l)/'2J 



,-,z) 




(14) 



Similarly, P{x, 0, z) can be expressed as 

- i^f}ir%:';:\y,)p(^'^i,o,.), (15) 

where f ^j is the binomial coefficient. Here, we also define (^j = 0, if x or y 
is not a non-negative integer, or if a; < y. Thus, the problem of calculating 
the first passage property is reduced to solving the self-consistency equations 
Eqs. (0) - (^). This is the key result in this section, and it will later 
serve as a basis for an iteration scheme. It is not unnatural that we end up 
with self-consistent relations rather than explicit solutions. The boundary 
conditions we have to satisfy at the interface between the two blocks are (1) 
continuity and (2) fiux conservation. Since these conditions are only implicit 
[i.e., they are relations among the fields P{x,y,z)], they result in implicit 
relations between G{x,a,z), which are the self-consistency conditions. 

Unfortunately, these conditions are essentially 4L — 2 coupled linear equa- 
tions, which are non-trivial to solve. We have developed an iterative scheme 
useful in getting the short time behaviors [n ~ L), which is discussed in 
Appendix A. We now develop an alternative method which can give the 
information about the asymptotic {n ^ L) behavior. 



3 Asymptotic Behavior 

3.1 Expansion of H+{z) 

We turn to an alternative method for obtaining the first passage time. We 
expand quantities in terms of the number of times the particle has crossed 
the interface between the blocks before reaching the boundary. Consider the 
+/— model again Eq. (Rf), and recall the previous definitions of P{x, y, z) and 



H^{z) (Eqs. (^) and (0)). We now define "half-space" Green's functions g'l 



as 



ry I / nn nn \ / /y nn' 



2 I \(x - a;' + y - y')/2/ Vl^^ " ^ + ^ + ^ )/2y 






2 [\{x' -x + y-y')/2j \{x' - x + y + y' - 2)/2J } 



where ("^j is defined as in Eq. (|T5|). Here, g'l {g'^) is the Green's function 
in the upper (lower) block with absorbing boundary a.t y = {y = 1). Due 
to the boundary condition, the functions (?^ do not contain contributions 
from the particles which leave the block, a property which will be useful 
subsequently. 

We define H^'{z), the part of H^{z) corresponding to particles which 
have crossed the interface n times. Using the definition of g^, 

Hf{z)=Y.zg^4L-l,y,Xo,yo,z), (17) 

y 

where we assume yo > without the loss of generality. We now calculate 
the flux of particles out of the upper block. Define P^"'\x,y,z) to be the 
part of P{x, y, z) corresponding to particles which have crossed the interface 
n times. At the edge of the upper block (the y = I line), 

P'''\xA,z)=g1{x,l,Xo,yo,z). (18) 

Half (1 — Py) of these particle will jump to (x + 1, 0). Therefore, the influx 
at the edge of the lower block {y = line), 

i^^)(a;,0,^) = ^pW(x-l,l,^). (19) 

We define the operator T±{x, x', z) as 

^-v{x ±l,z) = Y, r±(x, x\ z)v{x, /), (20) 

x' 

where v{x^ z) is a vector. Thus, Eq. (|TP|) in operator form is 

P{y,z)=T.{z)P^'\l,zl (21) 



(16) 



where we dropped the indices x and x'. 

We have to know what fraction the flux will go back to the upper block. 
We first obtain 



P^'\x,0,z) = Y.g'_{x,0,x',0,z)Pl^{x',0,z). (22) 

x' 

Again, half (py) of the particles will cross the interface and jump to (x— 1, 1), 
Pi^{x,l,z) ^ |p«(x+l,0,z) =T+(z)P«(0,z). (23) 

(2) 

Since H]_^ arises from from the walkers which have crossed the interface 

(2) 

twice, its sole contribution comes from -PjV. , which is 

H^^\z) = zJ2T.9liL-l,y,x',l,z)P^^ix',l,z). (24) 

y x' 

(2) 

We then calculate the fraction of P- which jumps back to the lower block, 
thus completing the cycle. At the edge of the upper block, 

P'^'\x,l,z) = Y.g'4x,l,x',l,z)Pl^{x,l,z), (25) 

x' 

and half of these will jump to {x + 1, 0) 

Pl^iO,z)=T.{z)P^'\l,z). (26) 

The above results can be written in a more compact form. We first define 
several operators 

{9l{z))x,x' = 9^{x,l,x',l,z), 
((?°(2;)),,,, = (?^(x,0,x',0,z), 

(hliz)),,,, = zJ2g^4L-l,y,x',l,z). (27) 

y 

in terms of which the above results can be written as 



pI^{1,z) = T^{z)g"_{z)T^{z)p('\l,z) 

in 



Hi'>{z) = h\{z)py^{l,z). (28) 

10 



Furthermore, by repeating the above procedure, we can show that 

Pl^\hz) = T^{z)g"_{z)T4z)gl{z)Pl^-'\l,z), 
H^:\z) = h\{z)P^f{l,z). (29) 



With one more definition 



u{z) = T+{z)g"_{z)T4z)gl{z), 
u^{z) = T+{z)g"_{z)T4z), (30) 

we arrive to the key result of this section: 

oo 

H4z) = H^'\z) +J2hl{zyiz)m{z)P^'\l,z). (31) 

The vahdity of the expansion has been checked by comparing H^{z) obtained 
above with that obtained by numerical simulations. Details of the simulations 
will be discussed later. 

3.2 Asymptotic form of H^ 

In this section, we derive the asymptotic form of the first passage time dis- 



tribution, starting from Eq. (|3TD. Recall the definition H^{z) = J^n.H'^z'^. 
In general, H^{z) is an infinite order polynomial of z, where if", the coeffi- 
cient of z'"-, is the hitting probability at time n. We now consider the various 
terms in Eq. (pi]). Using Eq. (|^), we can show the degree of H^'{z) can 
not be larger than 2L — 1. Since it does not give a contribution to H+{z) 
in the asymptotic n ^ L regime, we can ignore this term. Next, in the 
summand of the equation, the same operator u{z) is been repeatedly ap- 
plied to a vector u'{z)P^'^\l, z). Thus, we can approximate u^{z) with A*(z), 
where X{z) is the largest eigenvalue of u{z) If the operator in question is 
self-adjoint and diagonalizable, this approximation would surely be justified, 
at least for i 3> ll{\{z) — \2{z)), where \2{,z) is the second largest eigenvalue 
of u{z), but in this instance this step is an assumption, which however is 
supported by the numerical results below. We now ask whether the asymp- 
totic behavior will be changed by the approximation. The maximum degree 
of z for the terms in the summand can be calculated from Eq. (p!6[) . The 
maximum degree of h"^{z)^u{z)^u'{z) and P^'^\l^z) cannot be larger than 

11 



2L — 1, 4L — 2, 2L, 2L — 2, respectively. In the sum, the term containing u^{z) 
contributes for n < i{2L — 1) + 6L — 3. Therefore, the terms for which the 
eigenvalue approximation is not valid {i < io, and io is finite) give a contri- 
bution only up to finite time, and will not change the asymptotic behavior. 
Thus, 

oo 

HA^) - E ^^(^) ■ K{z)u'{z)P^'\l, z)]. (32) 

j=0 

The product h'^{z)u'{z)P^^^{l, z) is also a polynomial of finite order, and 
ignoring the product changes only the amplitude of the asymptotic behavior. 
We thus arrive to a simple expression 

-| oo 

^+W-r^T3 = E^^W- (33) 

We define the coefficient of the z^ term of X{z) to be cj. We can interpret 
X{z) as a generating function for a random walk process — the probability to 
jump j steps forward is given by the coefficient c]. The fraction of random 
walkers which survive after one step is So = A(l). The average displace- 
ment after one step is Si = A'(l)/A(l), and the average of the square of the 
displacement after one step is S2 = A"(l)/A(l), where 

A'(^) = z^Xiz), 
az 

\"{z) ^ [z^TKz). (34) 

We also define the variance o"^ = s^ — S2. Following the interpretation, the 
coefficient c* of the z^ term for A*(z) forms the distribution of the displace- 
ment of the random walker after i steps. The fraction of random walkers 
which survive for i steps is s^, the average displacement is isi, and variance 
ia"^. Since the second moment S2 is finite, we can apply the central limit 
theorem. Thus for large i, the coefficient c* becomes 

Substituting this to the equation for H^{z) Eq. (P^D, we obtain 

Hlr^y—^exp\- ^''''l'^' ], (36) 

12 



which can be evaluated by the method of the steepest descent to be 

ri n l-\nSo-{(T/sif/2 

H^ ~ exp Inso . , .^ 

si l-\nSo-{cr/siy 

= exp[— c(L)n], (37) 

where c{L) is size dependent decay constant. Even though the equation 
is derived for the +/— model, its derivation only assumes the existence of 
the half space Green's functions similar to Eq. (|16|). Since these functions 
exist for the most general situation of the two block system, the asymptotic 
distribution is a always simple exponential. Below, we compare the above 
results with an exact enumeration method, and find good agreement. 

3.3 Estimation of the Eigenvalue 

The problem of finding the asymptotic behavior of the first passage time 
distribution is reduced to finding X{z), the largest eigenvalue of the opera- 
tor u{z). Unfortunately, there is no known method to calculate the analytic 
expression of the eigenvalues of an arbitrary matrix, and the complicated 
structure of u{z) does not help matters. We present two methods to esti- 
mate A(z). These methods are not expected to produce exact numbers, but 
are intended to give some idea of the parameter L-dependences of the first 
passage time distribution. 

The first method is to express X{z) in terms of the average of the elements 
oiu{z). We start from the matrix T^[z)g'^{z), whose largest eigenvalue A+(z) 
is approximated as 

A+(^)~^7^ E E {T^{z)gl{z)).,.,, (38) 

where 2L — 1 is the size of the matrix. This approximation is motivated from 
the numerical fact that the eigenvector v^{z) corresponding to X+{z) is close 
to be uniform, i.e., {v+{z))x = Vo{z) for all x. Note that if the eigenvector is 
uniform, Eq. (|38|) becomes exact. We then obtain 

13 



The expression can be further siniphfied to 

A.W-| + ^|:^-*-"^(l4). (40) 

where we assume L ^ 1. The largest eigenvalue A_(z) of the matrix 
T_(z)(7° (z) can also be estimated by the same method. It turns out 

A_(^)~A+(z). (41) 

The matrix u{z) is given by the product of the two matrices, T+{z)g'L{z) and 
T_(z)(7" (z). We further approximate the largest eigenvalue of the product of 
two matrices as product of the largest eigenvalues of the two matrices, which 
implies 

Xiz) r^ X4z)X^iz) r^ Xliz). (42) 

We calculate A+(l), A'_|_(l) and A^(l), where the primed values are defined as 
the same way in Eq. (Q). By evaluating the integral in Eq. (0), we obtain 

A+(l) - (l: + ^]-^L-y' + 0{l/L), 



a;(i 







(43) 



The second approximation method is based on the interpretation that 
X{z) is related to a certain generating function for a random walk. Consider 
the matrix T_(l)(7° (1). The matrix gives the probability to reach points on 
y = 0, starting from points on y = 1 with an absorbing boundary at y = 0. 
Thus, X^{z) is roughly the generating function of the hitting probability on 
the line y = for a walk starting from y = 1. For simplicity, consider a walker 
starting from (0,1). Since the only effect of convection in the x direction 
is to remove all the walkers which do not reach y = until time step L, 
we only have to deal with an one-dimensional problem. The corresponding 
one-dimensional problem is treated by ignoring the x axis and limiting the 
maximum time step to be L. Then X-^-{z) is 

X+iz) ~ / rfn-==(l - exp-«/"), (44) 

Ji v27m 
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where the integrand is the flux to y = at step n, and we have approximated 
the sum by an integral. The eigenvalue A+(l) is 

/•oo 1 /-oo 1 

A+(l)~/ rfn-==(l-exp-«/")-/ dn-==(l-exp-«/"). (45) 

The first integral is the probability to hit y = during infinite period of time, 
which is unity. After simplifying the second integral for L ^ 1 we have 

The second approximation, compared to Eq. (|^), has the same dependence 
on L, but different numerical coefficients. This supports the previous sugges- 
tion that the coefficients obtained by these methods are not reliable. How- 
ever, the fact that two very different methods give the same dependence on 
L gives some support to the validity of the form. The leading term in A+(l) , 
which is the value of A+(l) in the limit L —>■ oo, deserves special attention. It 
is the probability that an unbiased random walker hits the y = line during 
an infinite period of time, which is equal to unity. Even in the case that 
the matrix is applied to the exact eigenvector, the random walker eventually 
has to be absorbed at the y = boundary for L — > oo, implying A+(l) = 1. 
Therefore, we set lim^^oo A+(l) = 1 from now on. 

We now proceed to calculate c(L). The expressions for So,Si and S2 can 
be obtained from Eqs. (|^) and (^21): 

So ~ A^(l)^l-^L-i/2 + 0(l/L), 



a' 




- -2^ + 2(^j ^j^L^^^^OiL), 



L"^ + 0{L). (47) 



150F 

Finally, we can calculate c{L) using the definition in Eq. 

In So 1 — lnSo(cr/si)^/2 



c{L) 



Si 1 -lnSo(or/Si)2 

-L-\ (48) 

32 ^ ' 
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3.4 Numerical Check 

We have calculated the decay constant c{L) in three major steps. First, 
we expand H^{z) in terms of the number of times a random walker have 
crossed the interface (Eq. (PT])). We then calculate its asymptotic distribution 
in terms of the eigenvalue X{z) with the help of the central limit theorem 
(Eq. (0)). We then estimate X{z) and c(L) (Eqs. (|7p - (|D). In this 
section, we check the validity of these results by comparing them with those 
of numerical simulations. It will serve as an intermediate check before we 
proceed to more general situation, where we will continue to use the methods 
developed above. 

We start with the exact sum Eq. (|31|). We calculate first few moments of 
H^{z) from the equation, which are easier to compare. The zero-th moment 
H+{1) is rather easy to evaluate. Consider the i-th term in the sum. We just 
have to calculate the product h'^{l)u{l)P^^ (1, 1), where we know all the 
individual terms. Furthermore, {i + l)-th term can be obtained by replacing 
pI^^^\i, 1) by i^^^*"*"^^(l, 1) = u{l)P^^^^\l, 1). We sum these terms in the 
increasing order of i, until the magnitude of the new term is smaller than 
certain value, which is chosen to be 10~^°. The higher moments are slightly 
more complicated to calculate. Consider again i-th term in the sum. The 
first moment can be calculated by using the chain rule 

[hl{z)u{z)Pl^^'\l,z)]' = hX{zMz)Pl^'''\l,z)+hl{z)u'{z)Pl^^'\l,z) 

+ hl{z)u{z)P!^'^'\l,z), (49) 



where the primed values are defined as the same way in Eq. (|34|) . Also 
{i + l)-th term can be obtained by replacing 



pS^'^(1,z) -. P^-\l,z)=u{z)P{^-''\l,zl 



m \ ' / m \ ' / \ ' m \ ; / ■ \ / m 



(50) 



Higher order moments are calculated by same way with the heavy use of the 
identity 

[A{z)B{z)p = ± (^^A^''\z)B^-'\z), (51) 
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where A^'^\z) is the n-th derivative of the function A(z). We compare the 



first five moments with those obtained by an exact enumeration method [|T5 
For several value of L = 10 ~ 100 and several initial conditions, the values 
obtained by the two methods are essentially identical. 

We check the next step, the asymptotic form of the hitting probability 
distribution Eq. (|37D. The form is simple exponential, and the decay constant 
c{L) is given in terms of X{z). We directly calculate the distribution iJ" by 
the exact enumeration method for L = 10 to 300. In Fig. 2, we show the 
distribution for L = 100. It is clear that iJ" is an exponential after some 
transition period. This is also true for the other sizes we have studied, and 
the length of the transition period is of the order L. 

We check the value of the decay constant c{L). Since the theoretical 
value of c{L) is given in terms of X{z), we have to know the value of X{z) 
in order to compare. To calculate the eigenvalue numerically, we go back 
to the discussion in the previous paragraph about calculating the moments 
of H+{z). We consider A(l) first. Since the matrix m(1) has been repeatly 



(^^+^)(1,1) to obtain i^^^+^)( 
A(l) = hm Pi^'-'\l, l)/Pi^^'\l, 1). (52) 



applied to the vector P- * (1, 1) to obtain P- * (1, 1 



We find the ratio hardly changes for i > 50, so we take the ratio at i = 100 
as A(l). For higher moments, start from the relation, 

X{z) = \imP^^^'\l,z)/P^^^'\l,z). (53) 

Taking the derivative and multiply z on both sides 

A'(l) = lim -^ ^ ' ^ m ^ ] [ m ^_^_m ^_2_ (54) 

(i^S (I'l))' 

The higher order terms (e.g., A"(l)) can be calculated in a similar way. Fi- 
nally, using these A(l), A'(l) and A"(l), the decay constant c(L) is calculated 
from Eq. (^). In Fig. 3, we show the values of c(L) just obtained as well as 
those obtained by numerical simulations, for several values of L. The simu- 
lational values are obtained by least square fitting the last 1/2 or 1/3 part 
of the numerically obtained if", like the one in Fig. 2. There are in general 
good agreements between these two values except for very small values of L, 
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where several approximations made to get the theoretical value may not be 
justified. 

We proceed to the last step of the calculation, the estimation of the decay 
constant c{L). In Fig. 4, we show the values of c{L) given by Eq. p8|) and 
those obtained by the exact enumeration. The values by the enumeration 
are identical to those shown in Fig. 3. It is clear that the numerical data 
shows the 1/L behavior as predicted by the theory. On the other hand, the 
measured prefactor (~ 1) is little smaller than the predicted value (39/32). 
These are all in accord to the expectation that the prediction of the 1/L 
dependence is reliable, but that of the prefactor is not. It is unexpected that 
the value of prefactor by the enumeration is so close to that of the theory. 

We have checked the steps to reach the decay constant c(L). The errors 
involved in the eigenvalue approximation are well controlled, and the approx- 
imation seems to be well justified for obtaining the asymptotic properties. 
Even though we do not have the same level of rigor in estimating the eigen- 
values X{z), we still have enough control to predict the correct dependence 
of L. 



4 The General Antisymmetric Model 

Having obtained a reasonable understanding of the asymptotic behavior of 
the +/— model, we turn to the general antisymmetric model. We now allow 
an arbitrary horizontal bias < p" < 1, and due to the symmetry in the 
system, we can restrict ourselves to p'^ > 1/2 without loss of generality. 

4.1 The Diffusive Limit: j9^ = 1/2 

We consider the antisymmetric model with no bias {p"^ = 1/2), first using the 
formalism developed for the +/— model. We now have a different operator 
u{z), and therefore a different eigenvalue X{z). The asymptotic distribution 
of i7^ is still a simple exponential, and the decay constant c{L,p'^) is given in 
terms of \{z) (Eq. (|48|)). In Sec. p.3| one estimate was based on by transform- 
ing the problem into an one dimensional diffusion problem with an absorbing 
boundary. In the transformation, one determines the average time required 
for the particles to be absorbed at the external boundaries at a; = ±L. In the 
+/— model, the transport in the horizontal direction is a purely convective 



process, so that the time is identical to the length L. Now, the transport in 

the horizontal direction is purely diffusive, so the time is 2L^. Substituting 

into Eq. (pi), 

2Q 
c{^p:)--^L-\ (55) 

The L~^ dependence is also consistent with the calculation for the one block 
case with no convection. In Fig. 5, we plot c{L,p^) determined by the above 
equation as well as those determined by the exact enumeration. The numer- 
ical data clearly shows the 1/L^ dependence with the prefactor about twice 
of that given in Eq. (|55D . 

A more direct check of these results is obtained by noting that in this 
case we are considering pure diffusion in a two dimensional strip of width 
2L, and one expects an exponential decay of tracer concentration with a time 
constant 0{L'^). More precisely, a straightforward solution of the diffusion 
equation for this geometry in Appendix B gives H^ ~ exp[— (7r^/8L^)?T,], in 
good agreement with the above simulation. 

4.2 The Neighborhood of j9" = 1 

It is useful to explicitly consider the case p" = 1 — e to first order of e and 
the formalism developed for the +/— model can be used with only minor 
changes. Starting from the master equation (|1|), it is straightforward to show 
that the half-space Green's functions g^,g'l Eq. (|^) must be modified to 

g1{x,y;x',y',z) = (1 - (x - x')e)(|r-^' 

J/ X — x' \ I X — x' \1 

[[{x-x' + y-y')/2) ~ [{x - x' + y + y')/2) j 

+ (x-x' + 2)e(^)^"^'+2 

' X- x' + 2 \ / X- x' + 2 \ 

^{x-x' + y-y' + 2)/2; ~ \{x - x' + y + y' + 2)/2j 

,z. 



X 



g1ix,y;x',y',z) = (1 - (x' - x)e){-r'' 



X 



x' — X \ I x' — X ^ 

(x' -x + y- y')/2 ~ [{x' - x + y + y' - 2)/2 , 
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U x'-x + 2 \ f x'-x + 2 ^ 

1 \{x' -x + y-y' + 2)/2 J Ux' -x + y + y')/2. 



(56) 



where terms of order e^ are ignored. The above equation reduces to Eq. (0) 
for e = 0, which suggests that the perturbation is not singular. The eigen- 
value X{z) for the matrix u{z) can be obtained by following the same pro- 
cedure as in the +/— model. We first calculate the eigenvalue \+{z) of the 
matrix T_{z)g'^{z) 

X+{z) ~ A+(2;)|,=o 



where \+{z)\^=q is the value of X+{z) at e = 0, and g°^{z) is defined as in 
Eq. (pTI). Replacing the sum by an integral, we obtain 

A+(l) ~ l-^(l-e)L-^/^ + C>(l/L) 



AV(1) - -^{l + e)L'/' + Oil) 

A'l(l) ~ j^{l + 3e)L'/' + 0{L'/% (58) 

where we have set lim^^^^oo A+(l) = 1 as discussed above. Similarly, the 
eigenvalue X^{z) of the matrix T+{z)g'L{z) is determined to be 

X4z)^X4z). (59) 

Combining these relations, the eigenvalue X{z) can be written as 

X{z) ~ X^{z)X+{z) ~ Xliz). (60) 

From X{z), the decay constant c{L,p^) is determined to be 

c{L,p:)r^-{l-2e)L-\ (61) 
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which is the same as the e = result aside from the factor 1 — 2e, whose origin 
is easy to understand. When we estimate the eigenvalue by mapping into an 
one dimensional problem, the average absorption time is required. For pure 
convection, this time is the length L divided by the horizontal velocity. In 
general, this velocity is given by 2p^ — 1 = 1 — 2e. Therefore, in the absence 
of diffusion, c(L,p") has to be modified by replacing L by L/{\ — 2e), which 
is precisely what we find by the expansion. In Fig. 6, we show c{L,p^) for 
L = 100 and several values of p^, where c{L,p^) is divided by the value at 
p" = 1. We also show the corresponding result by the expansion Eq. (|6T|) . 
There is a good agreement between the two, and the expansion is not only 
valid near p" = 1 but there is no systematic deviation down to p" = 0.65. 



4.3 Scaling Argument for General p 



The half space Green's functions Eq. ( |T^ ) have served as a starting point 
in the calculation of the decay constant at p" = 1 and its neighborhood. 
For other values of p", unfortunately, we are unable to find them in closed 
form. Although we can still show that the asymptotic distribution of if" is a 
simple exponential, without explicit knowledge of the Green's functions, we 
no longer can use the same method to estimate the eigenvalue and the decay 
constant. 

However, since we know the behavior for the two extreme cases of the 
antisymmetric model (p" = 1/2 and 1), we can try to bridge the gap by a 
simple scaling argument. We propose a scaling ansatz 

c{L,p-)=L-' f{^l (62) 

where f{x) is a scaling function which satisfies 

f{x) ~ I ^ l""^] (63) 

•^ ^ ^ [ X if X > 1. ^ ' 

where we have defined a crossover length L* = l/{2p^ — 1). (Recall that 
the lattice spacing has been set to one.) To verify that the ansatz is in 
accord with previous results, note first that for p!^ = 1, L* = 1, and since 
we are interested in L 3> 1, ( |63| ) gives c{L,p^) ~ L~^ as before. Next, for 
p" = 1/2, L* — >• oo and c{L,p]^) ~ L~^, again in agreement. Finally, for 
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p" = 1 — e, L* = 1/(1 — 2e), and for small e, L/ L* ^ 1, so the decay constant 
becomes (1 — 2e)/L, which is exactly Eq. (0). To verify the ansatz away 
from the limiting cases, we have used numerical simulation. In Fig. 7, we 
show the rescaled c(L,p") obtained by exact enumeration vs. the rescaled L, 
for several values of L = 10 ~ 200 and p" = 0.5 ~ 1.0. The data collapse 
into one scaling curve, which approaches a constant as a; — > and is linear 
for large x, precisely as expected from Eq. (|63D . Thus the scahng ansatz 
provides an excellent description of the general antisymmetric model. 

5 Conclusions 

We have studied the first passage time distribution iJ" of a two layer system 
of width L, and determined its asymptotic form to be an simple exponential 
decay in time. For the special case of an antisymmetric model, the decay 
constant is calculated using several techniques, and is found to cross over from 
the expected L~^ behavior in the pure-diffusion regime to an L^^ behavior 
at high velocities. 

The origin of the L^^ behavior in the convective regime is not intuitively 
obvious to us. It arises as the result of two contributions — one L~^'^ factor 
from In So term, and another L~^'^ factor from the 1/si term. As discussed in 
Sec. ^, So is roughly an eventual absorption probability of a one dimensional 
random walk, and Si is mean distance traveled before the absorption. This 
differs from a naive expectation that the L~^ behavior results from the mean 
distance Si behaving as L in the convective regime. 

Evidently, we have only considered the more tractable special cases in a 
two-layer system. It would be desirable to go beyond the antisymmetric limit 
of zero average velocity. In terms of the p" — p^ plane, the antisymmetric 
model corresponds to the line p" = 1 — p^, and the (elementary) one block 
case to the line p^ = pt- Due to the symmetries, the remaining region is 
bounded by the two cases with 1/2 < p^ < 1. Of course, one would like to 
consider convection in two directions as well multiple layers, but these rather 
more difficult problems must await further work. 



22 



Acknowledgements 

This research was supported in part by the Department of Energy under 
grant DE-FG02-93-ER14327. 



Appendix A — Iterative Method for Short Time 
Behavior 

We discuss an iterative scheme to obtain an approximate solution of the 
self-consistency equations Eqs. (0) - (|15D, based on interpreting them as a 
recursion relation. We input an trial solution of P{x, 0, z) and P{x, 1, z) to 
the equations, and obtain a (hopefully) improved approximation. In princi- 
ple, we repeat this procedure, until it converges to the correct solution. 
We start from a trial solution 

P(°)(x,0,2) = 0, 

P(°)(x, 1,^) = 0, (64) 

where the superscript indicates the number of iterations. For simplicity, 
we set Xo = and ?/o = 1. As shown in Eq. ([T3|) , the hitting probability 
H^ is related to G^ via H^{z) = zG^{L — 1,1,2;). Also, using Eq. (p^), 
G+{L — 1, 1, 2;) can be written as 

G4L- 1,1, z) = ^ E z'^-'^' {Fix' + 1,0, z)-P{x'- 1,1, z)) 

^ x'=-L+l 

+ z^-\ (65) 

Combining these relations, we obtain 

Hf\z) = z''. (66) 

This rather trivial result is due to the fact that the presence of the second 
block is ignored in the O*'^ order approximation. 

The next order values of P{x,0,z) and P{x,l,z) can be calculated by 
inserting the trial values into Eqs. (|l^ - (|l^) to obtain 

P^^\x,0,z) = 

P^'\xA,z) = (J)^fll (67) 



'2 \x/2j' 
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where we define r) = 0, ii x or y is not a non-negative integer, or if x < y. 



Similarly, we obtain 



1 ,^2 1 ^,/ X' 



Using Sterling's formula, and replacing the sum by an integral, 

H^^\z)c^z\1-Ml'/^). (69) 

The above result is unphysical, since H^ becomes negative for large L. The 
deficiency is due to the fact that we only include the fiux out of the upper 
block and not the fiux into the block, while both fiuxes are of the same order 
of magnitude. This problem will be resolved in the calculation at next order. 
The 2'^^^ order iterations of P{x, 0, z) and P{x, 1, z) are 

2...mii„.i,2 [(X' - x)/2P2> \{x'-l)/2) 

- |t(|)"'((.-V;'l)/2)(|)"-((;"lV2)''™' 
Again, using Sterling's formula and replacing the sum by an integral. 

P^^Hx,0.z) = - r^' dx'z^^'-^ ^ 



TT Jmax(x,i) ^i^x' -x){x' - 1) 



P^''\xXz) = z\\—-— r dx' , ^ (71) 

V TTX Tl Jl ^ (^x - X'){x' - I) 



r(2) 



and H\ becomes 



h'-'Hz) = z^ 



dx 



27r Jl \Jx — 1 
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Z^ /-L-l rx-l 1 

+ -— I dx I ax - 



27r^2 A /(a;_i_a;/)(a;'_ 1 



7^ /-O /-L-l , 1 

+ — dx dx'z'^-^^ 



27T J-L 



-1 Jl 



J{x' — a; — 1) (a; — 1) 



+ — / dx dx'z'^-'^ , (72) 

2vrA A+i J(^x' -x-l){x-l) 



Although the integrals in the equation can not be evaluated in closed form, 
we can understand their structure. The first term is the contribution of tracer 
which did not cross the boundary. The second and the third terms are amount 
of flux going out of the first block. Thus, these three terms are proportional 
to z^. The fourth and fifth terms are contributions from walkers which return 
to the upper block. Since the time to reach the boundary (x = L) depends 
on the where the walker exits {x') and reenters (x) the upper block, these 
terms contain different orders of z. 



Appendix B — The Single-layer System 

To provide some feeling for the more difficult case of a non-zero average 
velocity, we calculate the first passage time distribution for a single layer. 
Consider tracer moving between adsorbing boundaries at a; = ±L in the 
presence of a constant velocity field vx. The motion in the y-direction is 
simple diffusion, completely decoupled from that along x, and the problem 
effectively is one-dimensional. We have 

with the simple initial condition c{x, 0) = 6{x). Taking the Laplace transform 
via J^ dt e~^^, we have 

f)c 3 c 

.,.-Si.)^.-^D-, (74) 



which is readily solved for x 7^ as 

/ -T* 

(75) 



c^{x,s) = A^e^^/^sinh 
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p' + ailT^) 



Here the superscripts refer to x > and x < 0, respectively, and we have 
defined p = vL/2D and a = sL'^/D. The coefficients A^ are determined 
by the conditions c+ = c^ and dc^ /dx — dc~ /dx = —1/D at x = 0, which 
follow from the differential equation, so that 

A+ = A- = ^ , — ■ (76) 

ID-J-p^ + a cosh vp + (y 

The Laplace transform of the ffux leaving the system at x = ±L, which is 
identical to the ffrst passage time probability distribution, is 

J^is) = -D \ ' = , , , 77 

^ ' dx 2 cosh Vp^ + (T 

The long-time asymptotic behavior of J^ is controlled by the right-most 
singularities of the Laplace transform in the complex-s plane, in this case the 
poles where \/p^ + a = ±i7r/2 or s = s* = — 7r^D/4L^ — v'^ /AD. Thus, for 
t ~ oo, J^{t) ~ e''**. In the pure-diffusion limit, we set f = and recall that 
D = 1/2 and identify t with step number n, so that H/[ = J^{n) ~ e^'^ "/^■^ . 
In the opposite limit of large velocity, we see that J^{t) ~ e~^ t/io ^ which 
coincides with the long-time behavior at a fixed spatial point of the usual 
Gaussian solution of the CDE. 



26 



Figure Captions 

Fig. 1: System geometry: two semi-infinite blocks - y > and y < - with 
difi^erent velocities, with absorbing boundaries at x = ±L. 

Fig. 2: The hitting probability distribution H"^ for L = 100 obtained by the 
exact enumeration. An exponential behavior in the asymptotic regime 
is evident. 

Fig. 3: The decay constant c{L) for the +/— model given by the eigenvalue 
approximation (solid line) compared to the values obtained by the exact 
enumeration (diamonds). 

Fig. 4: The decay constants c{L) for the +/— model given by Eq. (^) 
(dashed line) and those obtained by the exact enumeration (diamonds). 

Fig. 5: The decay constants c{L,p'^) for the antisymmetric model with p^ = 
1/2 given by Eq. ( ^8]) (solid line) and those obtained by the exact 
enumeration (diamonds). 

Fig. 6: The normalized decay constants c(L,p") for the symmetric model 
predicted by the theory (solid line), and data obtained by the enumer- 
ation for L = 100 and several values of p" (diamonds). There is good 
agreement even down to p^ = 0.65. 

Fig. 7: The rescaled decay constants c{L,p^) obtained by the exact enumer- 
ation vs. rescaled L. Data for several different values of L = 10 ~ 200 
and p" = 0.5 ~ 1.0 collapse into one curve, which is the scaling function 

fix). 
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